#ifndef __fields
#define __fields
#include <iostream>
#include <vector>
#include "type.h"
#include <algorithm>
#include "stdio.h"
#include <fstream>
#include <sstream>
#include <cmath>
#include "physicalconstants.h"
#include "string.h"

///__Simbox中grid上的二维场群
class __Fieldgroup{
    public:
        /// field成员
        vector<vector<__Fields> > member;
        /// 每个field对应的grid，此次貌似没有用到，其实应该有关联性的，这个需要补充和修正
        vector<vector<__Vect2<double> > >grid;
        /// xygridmax
        __Vect2<long> xygridmax;
        /// fieldgroup的网格最小grid坐标
        __Vect2<long> xygridmin;
        /// xy 方向的空间步长
        __Vect2<double> dxy;
        /// 边界条件类型 1: damping  0: reflecting 位置：0, 下， 1,上， 2,左，3右        
        int bd_type[4];
        /// damping length 用于吸收边界的pml吸收厚度
        int damp_len[4];
        /// damping　开始的网格
        int startindex[4];

        __Fieldgroup(){}
        __Fieldgroup(const __Vect2<long>& xygridmax):xygridmax(xygridmax){
            xygridmin.member[0] = 0;
            xygridmin.member[1] = 0;
        }
        __Fieldgroup(const __Fieldgroup& other){
            this->xygridmax = other.xygridmax;
            this->xygridmin = other.xygridmin;
            this->grid = other.grid;
            this->member = other.member;
            this->dxy = other.dxy;
            memcpy(bd_type, other.bd_type, sizeof(int) * 4);
            memcpy(startindex, other.startindex, sizeof(long) * 4);
            memcpy(damp_len, other.damp_len, sizeof(int) * 4);
            
        }
        __Fieldgroup(const long& xgrid, const long& ygrid, __Vect2<double> dxy, const int bd_typer[]\
                , const int startindexr[], const int damp_lenr[]):dxy(dxy){
            //cout<<"grid max = "<<xgrid<<"\t"<<ygrid<<endl;
            xygridmax.member[0] = xgrid;
            xygridmax.member[1] = ygrid;
            xygridmin.member[0] = 0;
            xygridmin.member[1] = 0;
            memcpy(bd_type, bd_typer, sizeof(int) * 4);
            memcpy(startindex, startindexr, sizeof(int) * 4);
            memcpy(damp_len, damp_lenr, sizeof(int) * 4);
            //cout<<damp_len[0]<<"\t"<<damp_len[1]<<"\t"<<damp_len[2]<<endl;
            //cout<<damp_lenr[0]<<"\t"<<damp_lenr[1]<<"\t"<<damp_lenr[2]<<endl;
            //getchar();
        }
        ~__Fieldgroup(){}
        /// 初始化场群
        void init();
        /// 设置xygrid处的场
        void set_value(const __Vect2<long>& xygrid, __Fields& thisfield);
        friend ostream& operator<<(ostream& out, const __Fieldgroup& thisgroup){
            out<<thisgroup.xygridmin<<"\t"<<thisgroup.xygridmax<<"\t";
            for(int i = thisgroup.member.size(); i > 0; i--){
                for(int j = thisgroup.member[i].size(); j > 0; j--){
                    out<<thisgroup.member[i][j]<<"\t";
                }
                out<<endl;
                out<<endl;
            }
            return out;
        }
        /// 先damping B，半步推动B，再加damping E
        void advancepre(int rank, double delta_t);
        /// 推动一步E
        void advancenex(int rank, double delta_t);
        /// 输出数据
        void dump(const int& rank, int n, string fname, int sp);
        /// 设置 n类 xygrid处场数据 0 E, 1, B, 2 J
        void set_nvalue(int n, __Vect2<long>& xygrid, __Vect3<double>& myfield);
        /// Ampere定律更新电场
        void update_efield(double delta_t);
        /// Faraday定律更新磁场
        void update_bfield(double delta_t);
        /// 施加反射边界，如果bdtype已经设置，否则不起作用
        void apply_boundary();
        /// 场damp 乘， bd_type是1会起作用，type0：电场，1：磁场，axis：0 along x，1：along y， direction：0：从小到大，1：从大到小
        void field_damp_multipy(const int& bd_type, const int type, const int axis, const int direction, \
        const int& startindex, const int& damp_len);
        /// 场damp 除， bd_type是1会起作用，type0：电场，1：磁场，axis：0 along x，1：along y， direction：0：从小到大，1：从大到小
        void field_damp_divide(const int& bd_type, const int type, const int axis, const int direction, \
        const int& startindex, const int& damp_len);
        // 通过n设置特定的场，0--电场，1--磁场，2--电流密度
};

void __Fieldgroup::dump(const int& rank, int n, string fname, int sp){
    // ex
    // 记住此处c_str()函数用法 
    stringstream tempstream; 
    string tempint;
    if(n < 10){
        tempstream<<"0000";
    }
    else if(n < 100){
        tempstream<<"000";
    }else if(n < 1000){
        tempstream<<"00";
    }else if(n < 10000){
        tempstream<<"0";
    }
    tempstream<<n;
    tempstream<<"_rank_";
    tempstream<<rank;
    tempstream>>tempint;
    fname += tempint;
    /***
    tempstream.str("");
    // 先清空，后清空标志位
    tempstream.clear();
    tempint.clear();
    ***/
    int rankyj = rank / (rank / 2);
    int dumpstart = 2;
    int dumpend = 2;
    switch (sp)
    {
        case 0:
            {
                string filename = fname + "_rho.txt";
                ofstream myout;
                myout.open(filename.c_str());
                for(int i = member[0].size() - dumpend; i > dumpstart - 1; i -- )
                {
                //cout<<"dumping "<<filename<<endl;
                    for(int j = dumpstart; j < member.size() - dumpend + 1; j ++)
                    {
                        myout<<member[j][i].rho<<"\t";
                    }
                    myout<<endl;
                }
                myout.close();
                break; 
            }
        case 1:
        case 2:
        case 3:
            {
                string filename;
                if(sp == 1)
                {
                    filename = fname + "_ex.txt";
                }
                if(sp == 2)
                {
                    filename = fname + "_ey.txt";
                }
                if(sp == 3)
                {
                    filename = fname + "_ez.txt";
                }
                ofstream myout;
                myout.open(filename.c_str());
                //cout<<"dumping "<<filename<<endl;
                for(int i = member[0].size() - dumpend; i > dumpstart - 1; i --)
                {
                    for(int j = dumpstart; j < member.size() - dumpend + 1; j ++)
                    {
                        myout<<member[j][i].efield.member[sp - 1]<<"\t";
                    }
                    myout<<endl;
                }
                myout.close();
                break; 
            }
        case 4:
        case 5:
        case 6:
            {
                string filename;
                if(sp == 4)
                {
                    filename = fname + "_bx.txt";
                }
                if(sp == 5)
                {
                    filename = fname + "_by.txt";
                }
                if(sp == 6)
                {
                    filename = fname + "_bz.txt";
                }
                ofstream myout;
                myout.open(filename.c_str());
                //cout<<"dumping "<<filename<<endl;
                for(int i = member[0].size() - dumpend; i > dumpstart - 1; i --)
                {
                    for(int j = dumpstart; j < member.size() - dumpend + 1; j ++)
                    {
                        myout<<member[j][i].bfield.member[sp - 4]<<"\t";
                    }
                    myout<<endl;
                }
                myout.close();
                break; 
            }
        case 7:
        case 8:
        case 9:
            {
                string filename;
                if(sp == 7)
                {
                    filename = fname + "_jx.txt";
                }
                if(sp == 8)
                {
                    filename = fname + "_jy.txt";
                }
                if(sp == 9)
                {
                    filename = fname + "_jz.txt";
                }
                ofstream myout;
                myout.open(filename.c_str());
                //cout<<"dumping "<<filename<<endl;
                for(int i = member[0].size() - dumpend; i > dumpstart - 1; i --)
                {
                    for(int j = dumpstart; j < member.size() - dumpend + 1; j ++)
                    {
                        myout<<member[j][i].jfield.member[sp - 7]<<"\t";
                    }
                    myout<<endl;
                }
                myout.close();
                break; 
            } 
        default:
            break;
    }
}


void __Fieldgroup::set_value(const __Vect2<long>& xygrid, __Fields& thisfield){
    if (xygrid <= xygridmax && xygrid >= xygridmin){
        long i = xygrid.member[0], j = xygrid.member[1];
        this->member[i][j] = thisfield;
        //cout<<"pos now "<<this->member[i][j].position<<endl;
        this->member[i][j].position = xygrid;
        //cout<<"setting pos "<<this->member[i][j].position<<endl;
        //cout<<"will set field "<<i<<"\t"<<j<<"\t"<<thisfield<<endl;
        //cout<<"member now "<<member[i][j]<<endl;
        //getchar();
    }
}

void __Fieldgroup::set_nvalue(int n, __Vect2<long>& xygrid, __Vect3<double>& myfield){
    int i = xygrid.member[0], j = xygrid.member[1];
    //cout<<member[i][j].position<<"\t i, j "<<i<<j<<endl;
    //getchar();
    switch(n)
    {
        case 0:
            member[i][j].efield = myfield;
            break;
        case 1:
            member[i][j].bfield = myfield;
            break;
        case 2:
            member[i][j].jfield = myfield;
            break;
        default:
            break;
    }
}

void __Fieldgroup::init(){
    //cout<<"in init "<<xygridmax<<endl;
        //getchar();
    vector<__Fields> temp;
    //cout<<"initing "<<xygridmax<<endl;
    __Vect3<double> zeros(0, 0, 0);
    __Vect2<long> zeropos(0, 0);
    __Vect2<double> position(0, 0);
    __Fields voidfield(zeros, zeros, zeros, 0, zeropos);
    if(xygridmax >= xygridmin){
        for(int i = xygridmin.member[0]; i < xygridmax.member[0]; i++){
           vector<__Vect2<double> > thisx;
           position.member[0] = i * dxy.member[0];
           for(int j = xygridmin.member[1]; j < xygridmax.member[1]; j ++){
                position.member[1] = j * dxy.member[1];
                thisx.push_back(position);
                //cout<<"i, j"<<i<<"\t"<<j<<"\t"<<xygridmax<<endl;
                //getchar();
           }
           grid.push_back(thisx);
           thisx.clear();
        }
        // vector是pushback的先进去的用下标访问是对的
        // 我做过实验了。。。。
        for (int i = xygridmin.member[0]; i < xygridmax.member[0]; i++){
            //cout<<"i = "<<i<<endl;
            for(int j = xygridmin.member[1]; j < xygridmax.member[1]; j++){
                //cout<<"xygridmax:"<<xygridmax<<endl;
                //cout<<"j = "<<j<<endl;
                //cout<<"grid "<<grid[i][j]<<endl;
                voidfield.position.member[0] = i;
                voidfield.position.member[1] = j;
                temp.push_back(voidfield);
            }
            member.push_back(temp);
            temp.clear();
            vector<__Fields>(temp).swap(temp);
        }
    }
}

void __Fieldgroup::update_bfield(double delta_t){
    double dx = dxy.member[0];
    double dy = dxy.member[1]; 
    for(long i = 1; i< member.size() - 1; i++){
        //cout<<"i = "<<i<<endl;
        // 0 和最大grid处为guard，不需要计算
        for(long j = 1; j < member[i].size() - 1; j++){
             
            member[i][j].bfield.member[0] += - delta_t * \
                    ((member[i][j].efield.member[2] - member[i][j-1].efield.member[2])\
                    / dy);
            // dBx/dt = - dEz/dy
            // 这个地方如果i从1开始就会用到0的数据,就会出问题
            member[i][j].bfield.member[1] += delta_t * \
                    ((member[i][j].efield.member[2] - member[i-1][j].efield.member[2])\
                    / dx);
            // dBy/dt = dEz/dx
            member[i][j].bfield.member[2] += - delta_t * ((member[i+1][j].efield.member[1] - \
                        member[i][j].efield.member[1]) / dx - (member[i][j+1].efield.member[0]\
                            - member[i][j].efield.member[0]) / dy);
            /***
            member[i][j].bfield.member[0] += 0;
            // dBx/dt = - dEz/dy
            // 这个地方如果i从1开始就会用到0的数据,就会出问题
            member[i][j].bfield.member[1] += delta_t * \
                    ((member[i][j].efield.member[2] - member[i-1][j].efield.member[2])\
                    / dx);
            // dBy/dt = dEz/dx
            member[i][j].bfield.member[2] += - delta_t * (member[i+1][j].efield.member[1] - \
                        member[i][j].efield.member[1]) / dx;
            ***/
            
        }
    }
}

void __Fieldgroup::update_efield(double delta_t){
    double epsmu = 1.0 / (epsilon_0 * mu_0);
    double eps = 1.0 / epsilon_0;
    __Vect3<double> zeros(0, 0, 0);
    double dx = dxy.member[0], dy = dxy.member[1];
    for(long i = 1; i < member.size() - 1; i++){
        // 0 和最大grid处为guard，不需要计算
        for(long j = 1; j < member[0].size() - 1; j++){
            // 统一 012 xyz
            // 更新方法：半步B，一步E，再半步B
            // 这里使用了Ren.pdf上的算法，分TE模和TM模
            // dEx/dt = dBz/dy - Jx
            member[i][j].efield.member[0] += delta_t * ((member[i][j].bfield.member[2] - \
                        member[i][j-1].bfield.member[2]) * epsmu / dy - \
                    eps * member[i][j].jfield.member[0]);
            member[i][j].efield.member[1] += - delta_t * ((member[i][j].bfield.member[2] -\
                        member[i-1][j].bfield.member[2]) * epsmu / dx + \
                    eps * member[i][j].jfield.member[1]);
            
            // dEy/dt = - dBz/dx - Jy
            // 这里使用了Ren.pdf上的算法，分TE模和TM模
            member[i][j].efield.member[2] += delta_t * ((member[i+1][j].bfield.member[1] -\
                        member[i][j].bfield.member[1]) * epsmu / dx - \
                    (member[i][j+1].bfield.member[0] - member[i][j].bfield.member[0]) * \
                    epsmu / dy - eps * member[i][j].jfield.member[2]);
            /**
            //member[i][j].efield.member[0] += delta_t * ( - eps * member[i][j].jfield.member[0]);
            
            // dEy/dt = - dBz/dx - Jy
            member[i][j].efield.member[1] += - delta_t * ((member[i][j].bfield.member[2] -\
                        member[i-1][j].bfield.member[2]) * epsmu / dx + \
                    eps * member[i][j].jfield.member[1]);
            member[i][j].efield.member[2] += delta_t * (member[i+1][j].bfield.member[1] -\
                        member[i][j].bfield.member[1]) * epsmu / dx;
            ***/
        }
    }
}

void __Fieldgroup::advancepre(int rank, double delta_t){
    // 先damp b*
    /// 场damp 乘， bd_type是1会起作用，type0：电场，1：磁场，axis：0 along x，1：along y， direction：0：从小到大，1：从大到小
    
    field_damp_multipy(bd_type[0], 1, 0, 1, startindex[0], damp_len[0]);
    field_damp_multipy(bd_type[1], 1, 0, 0, startindex[1], damp_len[1]);
    field_damp_multipy(bd_type[2], 1, 1, 1, startindex[2], damp_len[2]);
    field_damp_multipy(bd_type[3], 1, 1, 0, startindex[3], damp_len[3]);
    ////cout<<rank<<" damping 3 "<<endl;
    //cout<<rank<<" b m finish "<<endl;
    // Faraday
    //cout<<"update b"<<endl;
    update_bfield(delta_t * 0.5);
       // 再 damp e* 
    //cout<<" damp e"<<endl;
    field_damp_multipy(bd_type[0], 0, 0, 1, startindex[0], damp_len[0]);
    field_damp_multipy(bd_type[1], 0, 0, 0, startindex[1], damp_len[1]);
    field_damp_multipy(bd_type[2], 0, 1, 1, startindex[2], damp_len[2]);
    field_damp_multipy(bd_type[3], 0, 1, 0, startindex[3], damp_len[3]);
    //cout<<rank<<" e m finish "<<endl;
	// Ampere
    //cout<<" update_efield "<<endl;
    //由此看来，对于MPI中，必须在计算完半步场后交换一次数据，才能继续计算，这场场才能传播过去
}

void __Fieldgroup::advancenex(int rank, double delta_t){

    // damp e/
    //cout<<rank<<" damp e"<<endl;
    field_damp_divide(bd_type[0], 0, 0, 1, startindex[0], damp_len[0]);
    field_damp_divide(bd_type[1], 0, 0, 0, startindex[1], damp_len[1]);
    field_damp_divide(bd_type[2], 0, 1, 1, startindex[2], damp_len[2]);
    field_damp_divide(bd_type[3], 0, 1, 0, startindex[3], damp_len[3]);
    //cout<<rank<<" e d finish "<<endl;
    // Faraday
    //cout<<" update_bfield"<<endl;
    update_bfield(delta_t * 0.5);
    // damp b/    
    //cout<<" damp b"<<endl;
    field_damp_divide(bd_type[0], 1, 0, 1, startindex[0], damp_len[0]);
    field_damp_divide(bd_type[1], 1, 0, 0, startindex[1], damp_len[1]);
    field_damp_divide(bd_type[2], 1, 1, 1, startindex[2], damp_len[2]);
    field_damp_divide(bd_type[3], 1, 1, 0, startindex[3], damp_len[3]);
    //cout<<rank<<" b d finish "<<endl;
    //cout<<" apply_boundary "<<endl;
    apply_boundary();
    //cout<<rank<<" field advance complete "<<endl;
    //cout<<"field advance complete "<<endl;
}

void __Fieldgroup::apply_boundary()
{
    // 0, 下， 1,上， 2,左，3右
    __Vect3<double> zeros(0, 0, 0);
    if(bd_type[0] == 0)
    {
        cout<<"down is reflecting"<<endl;
        for(int i = 1; i < member.size(); i ++)
        {
            member[i][1].efield = zeros;
        }
    }
    // 周期边界
    if(bd_type[0] == 2)
    {
        for(int i = 0; i < member.size(); i ++)
        {
            member[i][0] = member[i][member[0].size() - 2];
        }
    }
    if(bd_type[1] == 0)
    {
        cout<<"up is reflecting"<<endl;
        for(int i = 1; i < member.size(); i ++)
        {
            member[i][member[i].size()-1].efield = zeros;
        }
    }
    if(bd_type[1] == 2)
    {
        for(int i = 0; i < member.size(); i ++)
        {
            member[i][member[0].size() - 1] = member[i][1];
        }
    }

    if(bd_type[2] == 0)
    {
        cout<<"left is relfetcting"<<endl;
        for(int j = 1; j < member[0].size(); j ++)
        {
            member[1][j].efield = zeros;
        }
    }
    if(bd_type[3] == 0)
    {
        cout<<"right is reflecting"<<endl;
        for(int j = 1; j < member[0].size(); j++)
        {
            member[member.size() - 1][j].efield = zeros;
        }
    }
}



void __Fieldgroup::field_damp_multipy(const int& bd_type, const int type, const int axis, const int direction, \
        const int& startindex, const int& damp_len)
{
    __Vect3<double> zeros(0, 0, 0);
    //cout<<"bd_type = "<<bd_type<<endl;
    if(bd_type==1)
    {
    // direction 0: small to large, 1 large to small
    if(type == 0)
    {
        // efield
        // axis 0: boundary aligh x, 1: boundary align y
        for(int i = 1; i < damp_len; i ++)
        {
            int id = direction * (startindex - i) + (1 - direction) * (startindex + i);
            //cout<<id<<"\t"<<direction<<"\t"<<startindex<<endl;
            double damping = 1.0 - 0.5 * pow((1.0 * i) / damp_len, 3.0);
            /***
            cout<<" id "<<id<<endl;
            cout<<"damp_len "<<damp_len<<endl;
            cout<<"damping = "<<damping<<endl;
            ***/
            
            if(axis == 0)
            {
                for(int j = 1; j < member.size(); j ++)
                {
                    member[j][id].efield = member[j][id].efield * damping;
                }
            }
            else{
                for(int j = 1; j < member[1].size(); j ++)
                {
                    member[id][j].efield = member[id][j].efield * damping;
                }
            }
        }
    }
    else
    {
        for(int i = 1; i < damp_len; i ++)
        {
            int id = direction * (startindex - i) + (1 - direction) * (startindex + i);
            double damping = 1.0 - 0.5 * pow((1.0 * i) / damp_len, 3.0);
            /***
            cout<<" id "<<id<<endl;
            cout<<"damp_len "<<damp_len<<endl;
            cout<<"damping = "<<damping<<endl;
            cout<<" axis = "<<axis<<endl;
            ***/
            if(axis == 0)
            {
                for(int j = 1; j < member.size(); j ++)
                {
                    //cout<<"j id "<<j<<"\t"<<id<<endl;
                    member[j][id].bfield = member[j][id].bfield * damping;
                }
            }
            else{
                for(int j = 1; j < member[1].size(); j ++)
                {
                    //cout<<" id j"<<id<<"\t"<<j<<endl;
                    member[id][j].bfield = member[id][j].bfield * damping;
                }
            }
        }
    }
    }
}





void __Fieldgroup::field_damp_divide(const int& bd_type, const int type, const int axis, const int direction, \
        const int& startindex, const int& damp_len)
{
    __Vect3<double> zeros(0, 0, 0);
    if(bd_type==1)
    {
    // direction 0: small to large, 1 large to small
    if(type == 0)
    {
        // efield
        // axis 0: boundary aligh x, 1: boundary align y
        for(int i = 1; i < damp_len; i ++)
        {
            int id = direction * (startindex - i) + (1 - direction) * (startindex + i);
            double damping = 1.0 / (1.0 + 0.5 * pow((1.0 * i) / damp_len, 3.0));
            if(axis == 0)
            {
                for(int j = 1; j < member.size(); j ++)
                {
                    if(id > member[0].size())
                    {
                        cout<<" id over load "<<endl;
                    }
                    if(member[j][id].efield == zeros)
                    {}
                    else
                    {
                        member[j][id].efield = member[j][id].efield * damping;
                    }
                }
            }
            else{
                for(int j = 1; j < member[1].size(); j ++)
                {
                    
                    if(id > member.size())
                    {
                        cout<<" id over load "<<endl;
                    }
                    
                    
                    if(member[id][j].efield == zeros)
                    {}
                    else
                    {
                        member[id][j].efield = member[id][j].efield * damping;
                    }
                }
            }
        }
    }
    else
    {
        for(int i = 1; i < damp_len; i ++)
        {
            int id = direction * (startindex - i) + (1 - direction) * (startindex + i);
            double damping = 1.0 / (1.0 + 0.5 * pow((1.0 * i) / damp_len, 3.0));
            if(axis == 0)
            {
                if(id > member[0].size())
                {
                    cout<<" id over load "<<endl;
                }
                
                
                for(int j = 1; j < member.size(); j ++)
                {
                    if(member[j][id].bfield == zeros)
                    {}
                    else
                    {
                        member[j][id].bfield = member[j][id].bfield * damping;
                    }
                }
            }
            else{
                for(int j = 1; j < member[1].size(); j ++)
                {
                    
                    if(id > member.size())
                    {
                        cout<<" id over load "<<endl;
                    }
                    
                    if(member[id][j].bfield == zeros)
                    {}
                    else
                    {
                        member[id][j].bfield = member[id][j].bfield * damping;
                    }
                }
            }
        }
    }
    }
}



#endif
